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Abstract 

By using Moreau’s decomposition theorem for projecting onto cones, the 
problem of projecting onto a simplicial cone is reduced to finding the unique 
solution of a nonsmooth system of equations. It is shown that Picard’s 
method applied to the system of equations associated to the problem of pro¬ 
jecting onto a simplicial cone generates a sequence that converges linearly to 
the solution of the system. Numerical experiments are presented making the 
comparison between Picard’s and semi-smooth Newton’s methods to solve 
the nonsmooth system associated with the problem of projecting a point 
onto a simplicial cone. 
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1. Introduction 

The interest in the subject of projection arises in several situations, having 
a wide range of applications in pure and applied mathematics such as Convex 
Analysis (see e.g. ED, Optimization (see e.g. MUEmED, Numerical 
Linear Algebra (see e.g. 0), Statistics (see e.g. 0 EDI EH), Computer 
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Graphics (see e.g. [12] ) and Ordered Vector Spaces (see e.g. [131 EH 151 
EH EH EH]). More specifically, the projection onto a polyhedral cone, which 
has as a special case the projection onto a simplicial one, is a problem of 
high impact on scientific community^ The geometric nature of this problem 
makes it particularly interesting and important in many areas of science and 
technology such as Statistics (see e.g. HU), Computation (see e.g. inn, 
Optimization (see e.g. m e]) and Ordered Vector Spaces (see e.g. USD- 

The projection onto a general simplicial cone is difficult and computa¬ 
tionally expensive, this problem has been studied e.g. in [211 [22, '5][23,, 161I7]. 
It is a special convex quadratic program and its KKT optimality conditions 
form the linear complementarity problem (LCP) associated with it, see e.g 
PHEHH?]. Therefore, the problem of projecting onto simplicial cones can be 
solved by active set methods mmmm or any algorithms for solving 
LCPs, see e.g [25], [23] and special methods based on its geometry, see e.g 
[ 23! [23] . Other fashionable ways to solve this problem are based on the clas¬ 
sical von Neumann algorithm (see e.g. the Dykstra algorithm [231 [TUI 29]). 
Nevertheless, these methods are also quite expensive (see the numerical re¬ 
sults in {20J and the remark preceding section 6.3 in [3D]). 

In this paper we particularize the Moreau’s decomposition theorem for 
simplicial cones. This leads to an equivalence between the problem of pro¬ 
jecting a point onto a simplicial cone and one of finding the unique solution 
of a nonsmooth system of equations. We apply Picard’s method to find a 
unique solution of the obtained associated system. Under a mild assumption 
on the simplicial cone we show that the method generate a sequence that 
converges linearly to the solution of the associated system of equations. Nu¬ 
merical experiments are presented making the comparison between Picard’s 
and semi-smooth Newton’s methods for solving the nonsmooth system asso¬ 
ciated with the problem of projecting a point onto a simplicial cone. 

The organization of the paper is as follows. In Section [2] some notations, 
basic results used in the paper and the statement of the problems that we 
are interested are presented, in particular, the problem of projecting onto 
simplicial cone. In Section [3] we present some results about projection onto 
simplicial cones. In Section [4] we present two different Picard’s iterations for 
solving the problem of projecting onto simplicial cone. In Section [5] theoret- 


4 see the popularity of the Wikimization page Projection on Polyhedral Cone at 
www.convexoptimization.com/wikimization/index.php/SpeciahPopularpages 
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ical and numerical comparisons between Picard’s methods and semi-smooth 
Newton’s method for solving the problem of projecting onto simplicial cone 
[31J are provided. Some final remarks are made in Section [6] 

2. Preliminaries 

Consider R m endowed with an orthogonal coordinate system and let (■, •) 
be the canonical scalar product defined by it. Denote by || • || be the norm 
generated by (■,•). If a 6 K and x = (x 1 ,...,x m ) G M m , then denote 
a + := max{a,0}, a~ := max{—a, 0} and 

X+ := ((x 1 )^..., (x m ) + ) ,x~ := ((x 1 )-,..., , |x| := (|x J |,..., \x m \) . 

For x G M m , the vector sgn(x) will denote a vector with components equal to 
1, 0 or — 1 depending on whether the corresponding component of the vector 
x is positive, zero or negative. We will call a closed set K c M m a cone if 
the following conditions hold: 

(i) Ax + py G K for any A, p > 0 and x, y G K , 

(ii) x, —x G K implies x = 0. 

Let K c M m be a closed convex cone. The polar cone and the dual cone of 
K are, respectively, the sets 

K l = {x G | (x, y) <0, Vr/GiF}, K*:= {x G M m | (x, y) >0, Vr/G K}. (1) 

It is easy to see that /W = —K*. The set of all m x m real matrices 
is denoted by R mxm , I denotes the m x m identity matrix and diag(x) will 
denote a diagonal matrix corresponding to elements of x. 

For an M G R mxm consider the norm defined by 

\\M\\ := max{||Mx|| : x G M m , ||x|| = 1}, 
this definition implies 

||Mx||<||M||||x||, ||LM||<||L||||M||, (2) 

for any m x m matrices L and M. 
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Denote = (x = (x 1 ,..., x m ) E M m : X\ > 0,..., x m > 0} the nonneg¬ 
ative orthant. Let A E IR mxm be a nonsingular matrix. Then, the cone 

K:=AR™ = {Ax : x = (x\ ..., x m ) G M m , x x > 0,..., x m > 0}, (3) 

is called a simplicial cone or finitely generated cone. Let z G M m , then the 
projection Pk(z) of the point z onto the cone K is defined by 

P K (z) := argmin {\\z - y\\ : y G K} . 

From the definition of simplicial cone associated with the matrix A this 
definition is equivalent to 

Pk(z) argmin \\z — Ax || 2 : x — (x 1 ,..., x m ) GM m , Xi > 0,.. .,> x m > o|. 

Remark 1. It is easy to see that Pr™(x) = z + . It is well know that the 
projection onto a convex set is continuous and nonexpansive, in particular, 
we have ||z + — w + || < \\z — x|| for all z,w G M m , see fTj. 

The above remark shows that projection onto the nonnegative orthant is an 
easy problem. On the other hand, the projection onto a general simplicial 
cone is difficult and computationally expensive, this problem has been studied 
e.g. in [211 El [31, [7J [TBJ [22] . The statement of the problem that we are 
interested is: 

Problem 1 (projection onto a simplicial cone). Given A E M mxm a 
nonsingular matrix and z G R m , find the projection Pr(z ) of the point z 
onto the simplicial cone K = bUR’". 

The problem of projection onto a simplicial cone has many different formu¬ 
lations which allow us develop different techniques for solving them. In the 
next remark we present some of these formulations. 

Remark 2. Let A E M mxm be a nonsingular matrix and z G M m . From 
the definition of the simplicial cone associated with the matrix A in ([3]), 
the problem of projection onto a simplicial cone K = kUR™ may be stated 
equivalently as the following quadratic problem 

Minimize -\\z — Ax|| 2 , subject to x > 0. 
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Hence, ifv G R m is the unique solution of this problem then we have Pr(z) = 
v. The above problem is equivalent to the following nonnegative quadratic 
problem 


Minimize 


-x t Qx + x T b + c, 


subject to x > 0, 


(4) 


by taking Q = A T A, b = — A T z and c = z T z/2. The optimality condition 
for the problem ([4]) implies that its solution can be obtained by solving the 
following linear complementarity problem 


y = Qx + b, x > 0, y > 0, x T y = 0. (5) 


where y is a column vector of variables in M™. It is easy to establish that cor¬ 
responding to each nonnegative quadratic problems Q and each linear com¬ 
plementarity problems (J5]) associated to symmetric positive definite matrixes, 
there are equivalent problems of projection onto simplicial cones. Therefore, 
the problem of projecting onto simplicial cones can be solved by active set 
methods J25 . !ffl \2% \2ff or any algorithms for solving LCPs, see e.g [251 \2ffj 
and special methods based on its geometry, see e.g [25 l {&$ . Other fashionable 
ways to solve this problem are based on the classical von Neumann algorithm 
(see e.g. the Dykstra algorithm [28l I UTK 29]). Nevertheless, these methods 
are also quite expensive (see the numerical results in J20f and the remark 
preceding section 6.3 in 

As we will see in the next section, by using Moreau’s decomposition theo¬ 
rem for projecting onto cones, solving Problem [l] is reduced to solving the 
following problem. 

Problem 2 (nonsmooth equation). Given A G M mxm a nonsingular ma¬ 
trix and z G find the unique solution u of the nonsmooth equation 

(A T A — /) x + + x = A T z. (6) 

In this case, Pk{z) = Au + where K = AM™. 

Since x + = {x -\-\x\)/2 the Problem [2] is equivalent to the following problem: 

Problem 3 (absolute value equation). Given A G M mxm a nonsingular 
matrix and z G M m , find the unique solution u of the absolute value equation 

(A t A + I)x+(A t A-I)\x\=2A t z. (7) 

In this case, Pr{z) = Au + where K = AM™. 
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We will show in Section [4] that Problem [2] and Problem [3] can be solved by 
using Picard’s method. We end this section with the Banach’s hxed point 
theorem which will be used for proving our main result, its proof can be 
found in [32J (see Theorem 5.1 — 2 pag. 300 and Corollary 5.1 — 3 pag. 302). 


Theorem 1 (Banach’s fixed point theorem). Let(K,d ) be a non-empty 
complete metric space, 0 < a < 1 and T : X —» X a mapping satisfying 
d(T(x),T(y)) < ad(x, y ), for all x, yeX. Then there exists an unique x G X 
such that T(x) = x. Furthermore, x can be found as follows: start with an 
arbitrary element x 0 G X and define a sequence {x n } by x n+ \ = T(x n ), then 
lim„_ 5 . +00 x n = x and the following inequalities hold: 

cx. 

d(x,x n+ 1 ) < - d{x n+ i,x n ), d(x,x n+1 ) < ad(x,x n ), n = 0,1,.... 

1 — a 

3. Moreau’s decomposition theorem for simplicial cones 

In this section we present some results about projection onto simplicial 
cones. We recall the following result due to Moreau [33] : 

Theorem 2 (Moreau’s decomposition theorem). LetIi,L C M m be two 

mutually polar cones in M m . Then, the following statements are equivalent: 

(i) z = x + y, x € K, y € L and (x, y) = 0, 

(ii) x = Pk(z) and y = Pl{z)- 

Remark 3. Let K be a cone in W 71 . Note that from Moreau’s decomposi¬ 
tion theorem, definition of the polar cone and the dual cone in (JTJ) and the 
relationship iW = — K* it follows that 

P K {z) = z + P K * ( z)i V 2 e M m . 

Hence the problem of projecting onto K is equivalent to problem of projecting 
onto K*. 

The following result follows from the definition of the polar, see m- 
Lemma 1. Let A 6 ]g mxm a nonsingular matrix. Then, 

{AR ™) x = -(A t )- 1 M”\ 
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The following result has been proved in [13] by using Moreau’s decomposition 
theorem and Lemma Q] 

Lemma 2. Let A G 1Ri mxm be a nonsingular matrix and K = dR™ the 
corresponding simplicial cone. Then, for any z G W n there exists a unique 
x G such that the following two equivalent statements hold: 

(i) z = Ax + — (A T )- l x ~, x G M m , 

(ii) Ax + = Pk(z) and — {A T )~ l x~ = P K ±(z). 

The following result is a direct consequence of Lemma [2| it shows that solving 
Problem [l] is reduced to solving Problem [2] 

Lemma 3. Let A G M mxm be a nonsingular matrix, K = AM™ the corre¬ 
sponding simplicial cone and z G M m arbitrary. Then, equations i§ and 
0 have a unique solution u and P K (z ) = Au + , i.e., to solve Problem^is 
equivalent to solving either Problem or Problem [3| 

Proof. Since A is an m x m nonsingular matrix, multiplying by A T , the 
equality in item (i) of Lemma [2] is equivalently transformed into 

A T Ax + — x~ = A T z. 

As — x~ = x — x + , the above equality is equivalent to (pL Therefore, equation 
(|6]) is equivalent to the equation in item (i) of Lemma[2j Hence, we conclude 
from Lemma [2] that equation (6]) has a unique solution u and Pk{z ) = Au + . 
Since the equations © and (j7| are equivalent the result follows. □ 


4. Picard’s Method 

In this section we will present two different Picard’s iterations, one of 
them for solving Problem [2] and the other one for solving Problem [3} 

4-1. Picard’s Method for solving Problem [2] 

The Picard’s method for solving Problem [2] is formally defined by 

x k+ i = - (A T A - /) xf + A T z, k = 0,1,2,— (8) 

The sequence {£*,} with starting point x 0 G M m , called the Picard’s sequence 
for solving Problem [2j The next theorem provides a sufficient condition for 
the linear convergence of the Picard’s iteration (J8]) . 
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Theorem 3. Let A G ]^ mxm a nonsingular matrix, K = AMT the corre¬ 
sponding simplicial cone and z G M m arbitrary. If 

\\A r A-I\\<l, (9) 

then the Picard's sequence {x*,} for solving Problem^ converges to the unique 
solution u of equation (|6]) from any starting point xq G M m , Pr(z ) = Au + 
and the following error bound holds 

II U - x k \\ < ^ ^ ||x fc - Xfc-ill, V k — 1,2.... (10) 

Moreover, the sequence { x k } converges linearly to u as follows 

\\u - x k+ i\\ < \\A T A-I\\\\u - x k \\, k = 0,1,2,.... (11) 

Proof. Define the function F : M m —> M m as 


F(x) = — [A 1 A — /) x + + A T z. 


( 12 ) 


Since Remark [I] implies ||x + — y + 
easy to conclude that 


< ||x — y || for all x,y G M m , from (12) it 


\\F(x)-F(y)\\ < \\A T A — /||||x — y\\, Vx,?/GM m 

Therefore, as by assumption ||^4 T v4 — J|| < 1 we may apply Theorem [l] with 
X = M m , T = F, d(x,y ) = \\y — x|| for all x,y G M m and a = ||A T A — J||, for 
concluding that the Picard’s Method (|8]) or equivalently, the sequence 

x k+1 = F(x k ), k = 0,1,..., 


converges to a unique fixed point u of F, which from (12) is the solution of 
the Problem |2j i.e., 

(A T A - I)u + +u = A T z, 


and by using Lemma [3| w e have Pr{z ) = Au + . Moreover, Theorem [I] implies 
that the inequalities (10) and 0 hold. □ 






4-2. Picard’s Method for solving Problem [3] 

The Picard’s method for solving Problem [3] is formally defined by 

(A T A + l)x k+1 = - (A T A-I) \x k \ +2A r z, k = 0,1,2,.... (13) 


The sequence {a:*;} with starting point Xq G M m , called the Picard’s sequence 
for solving equation (J7]) or for projecting a point z G M m onto the simplicial 
cone K. From now on we will refer this method as Picard 2. 

Since A G M mxm is a nonsingular matrix we conclude that A T A is sym¬ 
metric and positive definite. Hence, A T A + / is nonsingular. Then for sim¬ 
plifying the notations define 

C := (A T A + I)- 1 (A T A - /). (14) 

Let Ai,... A m and ai,... a m be the eigenvalues of A T A and C, respectively. 
As A i > 0, for i — 1, 2,... m, it easy to conclude that 


Ill'll = max{|ay|,... \a m \} < 1, where 


<?i = 


1 - Aj 
Ai + 1 


i — 1, 2,... m. 


The next theorem provides the convergence of the Picard’s iteration (13). 


Theorem 4. Let A G M mxm be a nonsingular matrix, K = AM™ the corre¬ 
sponding simplicial cone and z G M m arbitrary. The Picard’s sequence 
for solving Problem [3] is well defined and converges to the unique solution 
u of equation ([7]) from any starting point x 0 G M m , Pr{z ) = Au + and the 
following error bound holds 


\u 


Xk\\ < 


lien 


lie|| 


||*^fc -^k— 1| 


V k = 1,2.... 


Moreover, the sequence {x k } converges linearly to u as follows 


(15) 


u-x k+1 \\ < ||e||||tt-a; fc n, k = 0,1,2,.... (16) 


Proof. Since the matrix A T A + I is nonsingular, the function F : M m —y M m , 


F(x) := - (A t A + /) 1 (A t A~ I)\x\ + 2(A t A + I) 1 A r z, (17) 


- 1 aT. 


is well defined. Since |||a;| 
(14) we conclude that 


< ||x — y || for all x,y G M m , from (17) and 


||F(a;)-F(y)|| < ||e|| \\x - y\\, Vx,ytR m . 
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Therefore, as ||C|| < 1 we may apply Theorem [I] with X = M™, T = F, 
d(x,y) = || y — x|| for all x, y G M™ and a = ||C||, for concluding that the 
Picard’s Method (13) or equivalently, the sequence 

x k+ i = F(x k ), k = 0,1,..., 


converges to a unique fixed point u of F, which from ( [T7| ) is the solution of 
the Problem [3j i.e., 

(. A T A + I)u+ (A T A - I) \u\ = 2A r z, 


and by using Lemma [3| w e have Pr{z) = Au + . Moreover, Theorem [I] implies 
that the inequalities (15) and (16) hold. □ 


5. Comparison between Picard’s and Newton’s methods 


In this section theoretical and numerical comparisons of above Picard’s 
methods and semi-smooth Newton’s method studied in 135 are provided. 


Also Picard’s method (13) is applied to solve an specific example. 


5.1. Theoretic comparison 

In this section theoretical comparisons between Picard’s methods and 
semi-smooth Newton’s method for solving Problem [l] will be provided. 

It is shown in [31] that the semi-smooth Newton method applied to equa¬ 
tion 0. namely, 

((A T A - /) diag(sgn(x£)) + I) x k+1 = A T z, k = 0,1,2,..., (18) 

is always well defined and under the assumption 

\\A t A - 7|| < b < t, (19) 

on the matrix A defining the simplicial cone K = AM™, the generated se¬ 
quence converges linearly to the unique solution u of Problem [2] from 
any starting point and, as a consequence of Lemma [3 we have Pr{z ) = Au + 
for any z G M m , which implies that u solves Problem l| 

Problem [lj i.e., the problem of projecting a point z G M m onto a simpli¬ 
cial cone K = AM™ is equivalent, by Lemma [3j to solving either Problem [2] 
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or 


Problem [3j Note that solving Problems [2] by Picard’s method (J8]) assump¬ 
tion ([9]) on the matrix A (see Theorem [ 3 ]) is less restrictive than assumption 
(19). When solving Problem [3] we only need the invertibility of the matrix A 
for Picard’s method ( fl3| to converge (see Theorem [4]). Therefore, Picard’s 
method (13) is theoretically more robust than Picard’s method ® and con¬ 


sequently than semi-smooth Newton method (18). In the next section we will 


present an example, where according to the established theory, only Picard’s 


method (13) can be applied. 


The main drawbacks of Picard (13) and semi-smooth Newton (18) is that 


both require the solution of a linear system in each iteration which constitute 
the largest computational effort of these methods. Picard’s method (|8]) do 
not have to solve a linear system, avoiding more complicated calculations, 
which is particularly interesting for large scale problems. We will investigate 


the efficiency of these methods in section 5.2.1 


5.1.1. Example 

Consider the monotone nonnegative cone, which is a simplicial cone K 
defined by 


K ■= {x = 


x 


,x m ) e 


X 


1 > X 2 > • • • > x m > 0} . 


( 20 ) 


The monotone nonnegative cone and the projection onto it occurs in various 
important practical problems such as the problem of map-making from rela¬ 
tive distance information e.g., stellar cartography (see web pag^Jand Section 
5.13.2 in [31j) and isotonic regression [35J |36] 3?i [38]- The isotonic regres¬ 
sion dH HH HD S23 is a very important topic in statistics with hundreds of 
papers and several books dedicated to this topic. This section provides a 
different view about projecting onto the monotone nonnegative cones via Pi¬ 


card’s method (13) which is related to the iterative theory of bidiagonal and 
tridiagonal matrices, and the Fibonacci numbers. The dual of the monotone 


5 www.convexoptimization.com/wikimization/index.php/Projection_on. 
Polyhedral_Convex_Cone 
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nonnegative cone is K* = AMT 1 , where A G M mxm is the nonsingular matrix 


( 1 \ 
-1 1 
-1 1 

V - 1 !/ 


A r A = 


/ 2 

-1 


-1 

2 

-1 


-1 


\ 


V 


2 -1 
-1 1 


From Remark |3j the problem of projecting a point onto K* is equivalent to 
projecting onto K. Let Ai,... X m be the eigenvalues of A T A. From [43J we 
have that the eigenvalues of matrix A T A are given by 


Ai 


2 + 2 cos 


2in 


2 m + 1 


‘8 = 1 , 2 ,..., m. 


( 21 ) 


Hence from (21) we conclude that 


0 < A* < 4, lim A m = 0, lim Ai = 4, 

m—>oo m—^oo 


lim ||A t A — ij| = 3. 

m—¥ oo 


Since ||A T A — J|| >1 for all m > 2, for projecting a point onto the cone 
K*, we can not apply semi-smooth Newton method studied in [31] neither 
Picard’s iteration (l8|). However Picard’s method (13) can be used. In order 


to reduce the computational cost of this method, for solving the linear system 
involved in each iteration, we suggest the following triangular decomposition 


A T A + / = 


/ d i -1 

G?2 


V 


\ 


-1 


dm— 1 


dr 


/ 




1 

d2 


\ 


-A l 


_ j_ 


V 


-A- 1 




where d m = 2, di = 3 — l/d i+ i for i = m — 1, • • • ,1. Another alternative for 
solving the linear system would be to compute the matrices 


R=(A t A + I) 1 (a t a~i), S=(A t A + I ) 1 a t . 


By using the recursion formulas for a tridiagonal matrix from [44], which are 
based on the results of [45; IS, HZ], after some algebraic manipulations and 
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taking into account that R is symmetric we obtain 


_ 2i 7 2iF 2m _2j+l 

F‘2rn +1 

F2tF‘2m.-2i ~ F2i-2F‘2m-2i+\ 


F, 


2m+l 


Rij = 


2-^2m-2j+l 

^2m+l 

2F 2l 

F 2 m+l 

2 

-^2m+l 

F2rn-2 


V F 


if 1 < i < j < m, 
if 1 < i = j < m, 
if 1 — i < j < m, 
if 1 < i < j = m, 
if i — 1, j — m, 
if i — j — 1, m. 


2m+l 


( F2iF 2m -2j+2 .r ■ 

— if * < J, 


s a = 


F 2m +i 

F 2 j-lF 2 m-2i+l 

F 2 m+1 

F 2 rn-2i+\ 

F 2m +1 


if 1 <j<i, 

if l — j <i. 


where F t is the Fibonacci sequence defined by F 0 — 0, F\ — 1 and F i+2 = 
Fi + F i+ 1- 

5.2. Computational results 

In this section we present two numerical experiments. In the first, nu¬ 
merical comparisons between Picard’s methods (J8]) , (13) and semi-smooth 
Newton’s method (18) for solving Problem [I] will be provided. In the sec¬ 
ond one, we study the behavior Picard’s method (13) solving the problem 
described in Section 5.1.1 All programs were implemented in MATLAB Ver¬ 
sion 7.11 64-bit and run on a 3A0GFfz Intel Core i5 — 4670 with 8.0 GB of 
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RAM. All MATLAB codes and generated data of this paper are available in 
http://orizon.mat.ufg.br/pages/34449-publications, 

General considerations: 

• In order to accurately measure the method’s runtime for a problem, 
each of them was solved 10 times and the runtime data collected. Then, 
we defined the corresponding method's runtime for a problem as the 
median of these measurements. 


• We consider that the method converged to the solution and stopped 
the execution when, for some k, the condition 


is satisfied. 



< RelativeTolerance , 


5.2.1. Numerical experiment I 

In this experiment, we study the percentage of problems for which a 
method was the fastest one (efficiency) to compare them. With the aim 
that methods (181), (13) and (18) find solutions on 1000 generated random 


test problems of dimension m = 1000, we construct the matrix A (defining 
the simplicial cone K = AM™) in each problem satisfying the condition ( [l9| ). 

We assume that a method is the fastest one for a problem, if the cor¬ 
responding runtime is less than or equal to 1.01 times the best time of all 
methods to find the solution. 

Each test problem was generated as follows: 


satisfying the condition (19), we 


To construct the matrix A G M. n ' 

Erst chose a random number b from the standard uniform distribution 
on the open interval (0,1/3). Then, we chose a random number b 
from the standard uniform distribution on the open interval (0, b). We 
computed the matrices S , V and D, respectively, from the singular value 
decomposition of a m x m generated real matrix containing random 
values drawn from the uniform distribution on the interval [—10 6 ,10 6 ]. 
Finally we computed 


A = S ( sqrt ( / + -V ) ) D, 


were v is the largest singular value of V and sqrt (I + ^V) is the square 
root of the diagonal matrix / + ^V. 
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We chose the solution u G M m containing random values drawn from the 
uniform distribution on the interval [—10 6 ,10 6 ] and computed 2 G 
from equation ([6]). Finally we chose a starting point Xq G M m containing 
random values drawn from the uniform distribution on the interval 

f-10 6 ,10 6 l. 


In order to provide information for the analysis of the large test problems 
set considered, we use the performance profiles (see [3Hj). The performance 
profile for a method is the cumulative distribution function for a performance 
metric. In this case we use the ratio of the method’s runtime versus the best 
runtime of all of the methods as the performance metric. Efficiency can be 
checked in the value of the profile function at 1. 

Figure [T| shows the performance profiles of the three methods for different 
relative tolerance values. These graphs reveal that Picard’s method (J8]) was 
the most efficient for low and medium accuracy, while semi-smooth Newton’s 


method (18) was the most efficient for high accuracy requirements. However, 


since semi-smooth Newton’s method (18) requires at each step the solution 


of a system of linear equations, which may become unreasonably expensive 
computationally as the problem dimension increases, these results suggest 
that for large scale problems Picard’s method pi) is recommended. 



(a) RelativeTolerance=10 7 


(b) RelativeTolerance=10 10 (c) RelativeTolerance=10 13 


Figure 1: Performance profiles on [1,4] for different accuracies. Picard, Picard2 and ss- 
Newton denotes the methods §), @ and (fl8|) , respectively. 


On the other hand, Picard’s method (13) was always the worst, except in 


the low accuracy case. It can be inferred from Figure [2j where convergence 
mean time for each problem consumed by Picard’s method (13) is less than 


consumed by semi-smooth Newton’s method (18). 
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Figure [2] shows, as one would expect, the number of iterations on semi¬ 
smooth Newton method is less than Picard’s methods (|8|) and (13) for solving 


the same set of problems and only for certain tolerance semi-smooth Newton 
method consumes less time. 



RelativeTolerance=10~ x , 7<x< 13 


RelativeTolerance=10" x , 7<x< 13 


Figure 2: Total overall iterations and total time in seconds, performed and consumed, re¬ 
spectively by each method to solve the 1000 test problems for different accuracies. Picard, 
Picard2 and ssNewton denotes the methods (|8|), (13) and (18), respectively. 


5.2.2. Numerical experiment II 

In this experiment, we study the behavior of Picard’s method (13) solving 
the problem described in Section 5.1.1 on sets of 100 generated random test 
problems of dimension m = 100, 500,1000,1500, 2000, respectively. 

Each m —dimensional test problem was generated as follows: We con¬ 
structed the matrix A (defining the simplicial cone K* = /UP™) as is defined 
in Section [5.1.1| We chose the solution u G M m , computed z G and chose 
a starting point xq G M m as we described in the previous Section |5.2.1 


The computational results obtained are reported in Table [TJ From these, 
it can be noted that for the same dimension, to achieve higher accuracy, the 
method needs to perform a greater number of iterations and consequently 
consume more runtime. The same behavior occurs when, for the same accu¬ 
racy, the dimension of the problem increases. 
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Dimension m 

Total Iterations 

Total Time 

100 

4927 

7475 

10036 

1.096 

1.624 

2.180 

500 

6613 

10333 

14055 

66.183 

103.411 

140.812 

1000 

8120 

12873 

17640 

449.507 

717.310 

984.274 

1500 

8159 

12924 

17732 

1358.698 

2151.743 

2952.247 

2000 

8814 

14054 

19359 

3098.215 

4955.041 

6820.121 

Relative Tolerance 

io- 7 

o 

1 

o 

o 

1 

CO 

o 

1 

-<1 

o 

1 

o 

o 

1 

CO 


Table 1: Total overall iterations and total time in seconds, performed and consumed, 
respectively by Picard’s method (13) to solve the 100 test problems of each dimension for 
different accuracies. 


6. Conclusions 


In this paper we studied the problem of projection onto a simplicial cone 
which, via Moreau’s decomposition theorem for projecting onto cones, is re¬ 
duced to finding the unique solution of a nonsmooth system of equations. 
Our main results show that, under a mild assumption on the simplicial cone, 
we can apply Picard’s method for finding a unique solution of the obtained 
associated system and that the generated sequence converges linearly to the 
solution for any starting point. Note that in Theorem [4] we do not make any 
assumption on the simplicial cone, on the other hand, we have to solve a 
linear equation in each iteration. It would be interesting to see whether the 
used technique can be applied for finding the projection onto more general 
cones. As has been shown in [7], the problem of projection onto a simplicial 
cone is reduced to a certain type of linear complementarity problem (LCP). 
Numerical comparisons between Picard’s methods ( 8fl3 ) and semi-smooth 
Newton’s method (18) for solving Problem [l] was provided in Section [5j It 
would also be interesting to compare these methods with the methods pro¬ 
posed in [221 EQj E] and the Lemke’s method for LCPs. 
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